Bayesian Regression Modeling to predict the third urine sample void volume and interval
Ted W. Simon
12/15/2020
Predicting Volume, Flow and Interval for the Third Sample
Here, Bayesian regression is used to predict the void volume, interval and urine production rate in the third sample from the 76 men in the NHANES sample from 2009-18 who provided three voids.
Also considered were the concentrations of sodium, potassium, chloride and bicarbonate in serum, anion gap, serum creatinine, and urinary creatinine. The use of these covariates,however, reduced the number of samples from 76 to 45 and for this reason, these were not used as covariates The models used are shown below.
First, though, some data wrangling is necessary.
men_all3 <- read_csv("men_all3_ions.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_double()
## )
## ℹ Use `spec()` for the full column specifications.dplyr::mutate(men_all3, grp = rep(0,76))## # A tibble: 76 x 22
## age wt ht bmi intrvl1 intrvl2 intrvl3 flow1 flow2 flow3 vol1 vol2
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 17 78.1 179. 24.3 150 75 56 0.007 0.4 3.36 1 30
## 2 11 28.3 142 14.0 160 37 23 0.138 0.243 0.739 22 9
## 3 6 22.2 120. 15.3 83 136 59 0.145 0.147 1.15 12 20
## 4 27 79.7 180. 24.6 119 40 60 0.319 0.15 0.283 38 6
## 5 7 34.7 127. 21.4 763 49 20 0.003 0.041 0.7 2 2
## 6 7 20.5 124. 13.4 51 78 36 0.49 0.09 0.639 25 7
## 7 8 48 139. 25.0 109 15 111 0.339 0.133 2.65 37 2
## 8 68 129. 182. 39.3 41 113 44 0.683 0.009 2.84 28 1
## 9 12 28 136. 15.2 41 19 59 0.463 0.579 2.37 19 11
## 10 10 38.6 135. 21.3 16 78 8 0.563 0.423 2.75 9 33
## # … with 66 more rows, and 10 more variables: vol3 <dbl>, S_osmol <dbl>,
## # SCr <dbl>, bicarb <dbl>, Cl <dbl>, Na <dbl>, Potass <dbl>, UCr <dbl>,
## # AUCr_pred <dbl>, grp <dbl>men_all3$grp <- rep(0,76)
for (i in 1:nrow(men_all3)) {
if (men_all3$intrvl1[i] <=120) { men_all3$grp[i] <- 1}
if (men_all3$intrvl1[i] >120 && men_all3$intrvl1[i] <=240) { men_all3$grp[i] <- 2 }
if (men_all3$intrvl1[i] >240 && men_all3$intrvl1[i] <=360) { men_all3$grp[i] <- 3 }
if (men_all3$intrvl1[i] >360 && men_all3$intrvl1[i] <=480) { men_all3$grp[i] <- 4 }
if (men_all3$intrvl1[i] >480) { men_all3$grp[i] <- 5 }
}
lmen3 <- data.frame(
log_age = log(men_all3$age),
log_wt = log(men_all3$wt),
log_ht = log(men_all3$ht),
log_bmi = log(men_all3$bmi),
log_int1 = log(men_all3$intrvl1),
log_int2 = log(men_all3$intrvl2),
log_int3 = log(men_all3$intrvl3),
log_vol1 = log(men_all3$vol1),
log_vol2 = log(men_all3$vol2),
log_vol3 = log(men_all3$vol3),
log_flow1 = log(men_all3$flow1),
log_flow2 = log(men_all3$flow2),
log_flow3 = log(men_all3$flow3),
grp = men_all3$grp )run models
if (file.exists("men_vol_mod_020321.Rds")) {
men_vol_model <- readRDS(file="men_vol_mod_020321.Rds")
} else {
men_vol_model <-
brm(data = lmen3, family = gaussian,
log_vol3 ~ (1 | grp) + log_age + log_wt + log_int1 + log_int2 + log_vol1 + log_vol2 + log_flow1 + log_flow2,
prior = c(set_prior("normal(0,1000)", class = "Intercept"),
set_prior("normal(0,100)", class = "b"),
set_prior("cauchy(0,2)", class = "sigma")),
warmup = 2000, iter = 5000, chains = 4,
control = list(adapt_delta = 0.99, max_treedepth=200))
saveRDS(men_vol_model, file="men_vol_mod_020321.Rds")
}
summary(men_vol_model)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: log_vol3 ~ (1 | grp) + log_age + log_wt + log_int1 + log_int2 + log_vol1 + log_vol2 + log_flow1 + log_flow2
## Data: lmen3 (Number of observations: 76)
## Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~grp (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.48 0.52 0.01 1.83 1.00 3708 5369
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 5.92 1.70 2.42 9.18 1.00 9579 8662
## log_age -0.46 0.27 -1.00 0.07 1.00 7912 7955
## log_wt 0.16 0.39 -0.60 0.92 1.00 8266 8210
## log_int1 -9.09 7.93 -24.62 6.51 1.00 4949 6705
## log_int2 45.12 19.33 6.24 82.26 1.00 5192 6886
## log_vol1 9.05 7.91 -6.46 24.48 1.00 4949 6704
## log_vol2 -45.42 19.30 -82.55 -6.61 1.00 5189 6814
## log_flow1 -9.14 7.99 -24.76 6.52 1.00 4946 6718
## log_flow2 45.56 19.36 6.57 82.77 1.00 5192 6776
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.97 0.09 0.82 1.15 1.00 10958 7967
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).if (file.exists("men_intrvl_mod_020321.Rds")) {
men_intrvl_model <- readRDS(file="men_intrvl_mod_020321.Rds")
} else {
men_intrvl_model <-
brm(data = lmen3, family = gaussian,
log_int3 ~ (1 | grp) + log_age + log_wt + log_ht + log_bmi + log_int1 + log_int2 + log_vol1 + log_vol2 + log_flow1 + log_flow2,
prior = c(set_prior("normal(0,1000)", class = "Intercept"),
set_prior("normal(0,100)", class = "b"),
set_prior("cauchy(0,2)", class = "sigma")),
warmup = 2000, iter = 5000, chains = 4,
control = list(adapt_delta = 0.99, max_treedepth=200))
saveRDS(men_intrvl_model, file="men_intrvl_mod_020321.Rds")
}
summary(men_intrvl_model)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: log_int3 ~ (1 | grp) + log_age + log_wt + log_ht + log_bmi + log_int1 + log_int2 + log_vol1 + log_vol2 + log_flow1 + log_flow2
## Data: lmen3 (Number of observations: 76)
## Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~grp (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.46 0.48 0.01 1.71 1.00 2439 4498
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 239.87 291.91 -335.45 810.02 1.00 4955 6338
## log_age 0.06 0.15 -0.23 0.36 1.00 12025 8305
## log_wt 25.95 31.73 -36.60 87.80 1.00 4959 6385
## log_ht -51.11 63.42 -174.92 73.83 1.00 4953 6362
## log_bmi -26.34 31.74 -88.35 36.33 1.00 4954 6411
## log_int1 -7.45 4.29 -15.75 1.15 1.00 5695 6472
## log_int2 5.42 10.70 -16.13 26.29 1.00 6313 6683
## log_vol1 7.58 4.28 -0.96 15.81 1.00 5680 6435
## log_vol2 -5.63 10.68 -26.48 15.94 1.00 6312 6590
## log_flow1 -7.67 4.32 -15.97 0.96 1.00 5681 6484
## log_flow2 5.59 10.72 -16.01 26.50 1.00 6314 6622
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.52 0.05 0.43 0.62 1.00 9371 8105
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).if (file.exists("men_flow-mod_020321.Rds")) {
men_flow_model <- readRDS(file="men_flow-mod_020321.Rds")
} else {
men_flow_model <-
brm(data = lmen3, family = gaussian,
log_flow3 ~ (1 | grp) + log_age + log_wt + log_ht + log_bmi + log_int1 + log_int2 + log_vol1 + log_vol2 + log_flow1 + log_flow2,
prior = c(set_prior("normal(0,1000)", class = "Intercept"),
set_prior("normal(0,100)", class = "b"),
set_prior("cauchy(0,2)", class = "sigma")),
warmup = 2000, iter = 5000, chains = 4,
control = list(adapt_delta = 0.99, max_treedepth=200))
saveRDS(men_flow_model, file="men_flow-mod_020321.Rds")
}
summary(men_flow_model)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: log_flow3 ~ (1 | grp) + log_age + log_wt + log_ht + log_bmi + log_int1 + log_int2 + log_vol1 + log_vol2 + log_flow1 + log_flow2
## Data: lmen3 (Number of observations: 76)
## Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~grp (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.44 0.48 0.01 1.73 1.00 3727 5742
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -101.68 341.17 -766.70 575.87 1.00 5979 6971
## log_age -0.51 0.28 -1.07 0.04 1.00 12271 7932
## log_wt -11.57 37.08 -83.84 61.97 1.00 5974 6886
## log_ht 22.22 74.11 -124.56 166.59 1.00 5979 6977
## log_bmi 12.33 37.09 -61.06 84.43 1.00 5974 6876
## log_int1 -1.92 8.18 -17.61 14.44 1.00 5479 6331
## log_int2 39.09 19.57 1.04 78.18 1.00 5237 6657
## log_vol1 1.78 8.16 -14.57 17.54 1.00 5480 6393
## log_vol2 -39.17 19.53 -78.17 -1.30 1.00 5239 6670
## log_flow1 -1.80 8.25 -17.63 14.67 1.00 5477 6290
## log_flow2 39.35 19.59 1.31 78.45 1.00 5240 6657
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.99 0.09 0.83 1.18 1.00 11777 8388
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).The desired result of these three models is the daytime within-person distributional parameters forf void volume, void interval and urine production rate. Each of these was modeled as a lognormal distribution with parameters mu and sigma, the log mean and log standard deviation respectively. For each individual, these two parameter are calculated from three measurements. The next chunk provides the calculations for these, both from the data and for the model predictions. Note that only the natural logarithms of the third void volume, void interval and urine production rate are the direct model outputs.
mutate(lmen3, mu_vol=rep(0,76), sig_vol=rep(0,76), mu_intrvl=rep(0,76), sig_intrvl=rep(0,76))## log_age log_wt log_ht log_bmi log_int1 log_int2 log_int3 log_vol1
## 1 2.833213 4.357990 5.189060 3.190065 5.010635 4.317488 4.025352 0.0000000
## 2 2.397895 3.342862 4.955827 2.641198 5.075174 3.610918 3.135494 3.0910425
## 3 1.791759 3.100092 4.790820 2.728506 4.418841 4.912655 4.077537 2.4849066
## 4 3.295837 4.378270 5.193512 3.201526 4.779123 3.688879 4.094345 3.6375862
## 5 1.945910 3.546740 4.845761 3.065725 6.637258 3.891820 2.995732 0.6931472
## 6 1.945910 3.020425 4.819475 2.591516 3.931826 4.356709 3.583519 3.2188758
## 7 2.079442 3.871201 4.932313 3.216874 4.691348 2.708050 4.709530 3.6109179
## 8 4.219508 4.862908 5.201256 3.670715 3.713572 4.727388 3.784190 3.3322045
## 9 2.484907 3.332205 4.910447 2.721953 3.713572 2.944439 4.077537 2.9444390
## 10 2.302585 3.653252 4.903050 3.057298 2.772589 4.356709 2.079442 2.1972246
## 11 3.496508 4.117410 5.178971 2.969902 5.111988 4.077537 4.691348 2.9957323
## 12 2.833213 4.117410 5.150397 3.026746 6.492240 3.555348 4.442651 2.7725887
## 13 3.496508 4.617099 5.234845 3.357594 4.158883 5.111988 3.135494 3.5263605
## 14 2.833213 4.339902 5.180097 3.190065 4.624973 3.737670 4.248495 3.5263605
## 15 2.397895 3.468856 4.971894 2.735665 5.062595 2.772589 4.219508 3.0910425
## 16 4.382027 4.453184 5.156178 3.349904 4.955827 4.718499 4.127134 2.7725887
## 17 1.945910 3.269569 4.810557 2.856470 3.583519 3.091042 4.465908 2.5649494
## 18 2.079442 3.169686 4.829912 2.721295 4.127134 4.382027 4.219508 2.4849066
## 19 4.025352 4.576771 5.170484 3.446808 4.905275 4.762174 4.356709 3.0910425
## 20 3.637586 4.363099 5.182345 3.206803 4.394449 3.713572 4.615121 0.0000000
## 21 2.197225 3.538057 4.936630 2.873565 6.278521 4.770685 3.713572 3.0445224
## 22 2.772589 4.099332 5.115596 3.077312 5.141664 3.871201 4.584967 3.7135721
## 23 4.343805 4.442651 5.182345 3.288402 4.700480 4.442651 3.850148 2.8903718
## 24 3.610918 4.290459 5.135798 3.230804 4.682131 4.204693 4.477337 3.0910425
## 25 2.944439 4.203199 5.227358 2.960105 3.713572 4.709530 3.850148 2.3025851
## 26 2.890372 4.007333 5.175019 2.867899 4.672829 3.850148 4.406719 3.3672958
## 27 4.382027 3.843744 5.148657 2.753661 4.158883 4.962845 3.526361 2.6390573
## 28 2.708050 3.749504 5.057519 2.844909 4.605170 4.859812 4.465908 2.7725887
## 29 1.945910 3.186353 4.799914 2.797281 4.262680 3.401197 3.688879 2.7080502
## 30 3.258097 4.269697 5.132853 3.214868 4.477337 4.955827 3.637586 3.7376696
## 31 3.637586 4.694096 5.246498 3.411148 3.610918 4.465908 3.891820 2.4849066
## 32 1.945910 3.238678 4.847332 2.753661 4.905275 3.555348 3.555348 3.2188758
## 33 4.025352 4.371976 5.139322 3.303217 5.236442 3.258097 4.356709 3.7612001
## 34 2.197225 3.198673 4.821893 2.766319 5.068904 3.610918 4.043051 3.4657359
## 35 2.079442 3.613617 4.906015 3.010621 3.931826 4.488636 3.871201 2.7080502
## 36 3.637586 4.278054 5.218191 3.054001 4.330733 4.477337 4.343805 2.8903718
## 37 3.931826 4.080922 5.153869 2.985682 3.401197 4.927254 3.737670 2.4849066
## 38 2.079442 3.360375 4.906755 2.760010 4.553877 4.356709 3.526361 2.7725887
## 39 2.484907 3.616309 4.980863 2.862201 3.761200 4.406719 3.401197 2.4849066
## 40 4.382027 4.163560 5.172754 3.030134 4.912655 4.828314 3.871201 1.7917595
## 41 2.079442 3.546740 4.862908 3.030134 4.890349 4.077537 4.234107 2.7725887
## 42 2.197225 3.152736 4.850467 2.660260 3.496508 3.465736 3.737670 2.3025851
## 43 3.465736 4.453184 5.177279 3.310543 5.003946 4.394449 2.564949 3.7135721
## 44 2.772589 4.175925 5.170484 3.044522 5.707110 3.496508 4.043051 3.4011974
## 45 1.945910 3.095578 4.794964 2.714695 4.077537 3.663562 3.828641 2.6390573
## 46 4.174387 4.280824 5.139322 3.210844 4.852030 4.762174 4.094345 3.9889840
## 47 1.791759 3.144152 4.771532 2.809403 4.406719 3.806662 3.828641 3.1780538
## 48 1.791759 3.054001 4.805659 2.653242 4.158883 3.931826 4.605170 2.9444390
## 49 2.302585 3.481240 4.955827 2.778819 4.584967 3.332205 3.871201 3.8066625
## 50 2.890372 4.318821 5.144583 3.238678 3.912023 4.043051 4.532599 2.9444390
## 51 2.708050 4.261270 5.174453 3.122365 4.234107 4.574711 4.219508 3.3322045
## 52 2.197225 3.325036 4.964940 2.602690 3.970292 4.248495 3.951244 2.8903718
## 53 3.583519 4.701389 5.167639 3.575151 4.744932 3.663562 4.499810 3.4657359
## 54 3.465736 4.784153 5.182907 3.629660 4.532599 3.737670 2.564949 2.5649494
## 55 3.828641 4.235555 5.144000 3.157000 3.737670 4.330733 4.532599 1.9459101
## 56 3.871201 4.644391 5.207298 3.440418 4.691348 4.624973 3.663562 3.4339872
## 57 2.833213 4.107590 5.179534 2.960105 4.442651 3.951244 4.343805 2.8903718
## 58 2.397895 4.003690 5.049215 3.113515 4.828314 4.499810 3.295837 2.5649494
## 59 2.995732 4.251348 5.118592 3.222868 3.850148 4.382027 4.025352 2.3978953
## 60 1.791759 3.000720 4.736198 2.740840 4.262680 4.454347 3.178054 2.8903718
## 61 3.912023 4.077537 5.195731 2.895912 4.043051 4.605170 3.784190 2.7080502
## 62 2.397895 3.580737 4.909709 2.970414 4.919981 4.110874 4.025352 2.3978953
## 63 2.639057 4.200205 5.131672 3.148453 5.017280 0.000000 3.970292 3.6635616
## 64 2.564949 4.102643 5.118592 3.077312 3.988984 3.988984 4.304065 3.2958369
## 65 2.639057 4.135167 5.053695 3.238678 4.553877 4.158883 3.891820 3.1354942
## 66 1.945910 3.161247 4.849684 2.674149 4.663439 3.828641 3.988984 3.2958369
## 67 4.382027 4.347694 5.143416 3.269569 5.153292 4.442651 4.727388 3.4657359
## 68 4.077537 4.056989 5.083266 3.100092 4.820282 4.521789 4.718499 2.6390573
## 69 2.890372 4.070735 5.156178 2.970414 4.605170 4.025352 3.988984 2.3025851
## 70 2.772589 4.154185 5.171620 3.020425 4.595120 3.663562 4.634729 3.2580965
## 71 2.397895 3.575151 5.000585 2.785011 3.891820 4.574711 3.433987 1.9459101
## 72 4.110874 4.764735 5.203457 3.569533 5.789960 4.532599 3.663562 2.7080502
## 73 4.189655 4.453184 5.200153 3.261935 4.343805 4.804021 3.526361 3.1354942
## 74 2.484907 3.380995 4.950885 2.687847 3.044522 4.262680 4.219508 1.7917595
## 75 3.737670 4.309456 5.174453 3.169686 6.315358 4.882802 3.433987 3.4011974
## 76 3.761200 4.768988 5.202907 3.572346 4.043051 4.369448 4.343805 2.0794415
## log_vol2 log_vol3 log_flow1 log_flow2 log_flow3 grp mu_vol sig_vol
## 1 3.4011974 5.236442 -4.9618451 -0.91629073 1.21104772 2 0 0
## 2 2.1972246 2.833213 -1.9805016 -1.41469384 -0.30245736 2 0 0
## 3 2.9957323 4.219508 -1.9310215 -1.91732269 0.14236724 1 0 0
## 4 1.7917595 2.833213 -1.1425642 -1.89711998 -1.26230838 1 0 0
## 5 0.6931472 2.639057 -5.8091430 -3.19418321 -0.35667494 5 0 0
## 6 1.9459101 3.135494 -0.7133499 -2.40794561 -0.44785082 1 0 0
## 7 0.6931472 5.683580 -1.0817552 -2.01740615 0.97418221 1 0 0
## 8 0.0000000 4.828314 -0.3812604 -4.71053070 1.04415610 1 0 0
## 9 2.3978953 4.941642 -0.7700282 -0.54645280 0.86415498 1 0 0
## 10 3.4965076 3.091042 -0.5744757 -0.86038310 1.01160091 1 0 0
## 11 2.4849066 3.401197 -2.1202635 -1.59454930 -1.29098418 2 0 0
## 12 2.9444390 4.553877 -3.7297014 -0.61064596 0.11154137 5 0 0
## 13 2.1972246 3.332205 -0.6329933 -2.91877123 0.19638881 1 0 0
## 14 2.7080502 3.784190 -1.0996128 -1.03001950 -0.46362402 1 0 0
## 15 0.6931472 3.850148 -1.9732813 -2.07944154 -0.36961546 2 0 0
## 16 3.3672958 2.944439 -2.1803675 -1.35092722 -1.18417018 2 0 0
## 17 1.7917595 3.871201 -1.0188773 -1.29828348 -0.59420723 1 0 0
## 18 3.2958369 4.564348 -1.6398971 -1.08470938 0.34500714 1 0 0
## 19 2.1972246 1.386294 -1.8140051 -2.56394986 -2.97592965 2 0 0
## 20 0.6931472 5.003946 -4.4228486 -3.01593498 0.38865799 1 0 0
## 21 3.2188758 4.442651 -3.2441936 -1.55116900 0.72899683 5 0 0
## 22 0.0000000 4.574711 -1.4271164 -3.86323284 -0.01005034 2 0 0
## 23 2.4849066 1.945910 -1.8078889 -1.95899539 -1.90380897 1 0 0
## 24 3.0445224 4.007333 -1.5896353 -1.16155209 -0.47000363 1 0 0
## 25 3.1354942 2.302585 -1.4105871 -1.57503649 -1.54646311 1 0 0
## 26 2.7725887 3.135494 -1.3056365 -1.07880966 -1.27296568 1 0 0
## 27 3.2958369 2.484907 -1.5186835 -1.66600826 -1.04128722 1 0 0
## 28 3.7376696 3.663562 -1.8325815 -1.12085790 -0.80296205 1 0 0
## 29 3.0445224 5.036953 -1.5558971 -0.35667494 1.34807315 1 0 0
## 30 2.7725887 1.609438 -0.7402388 -2.18036746 -2.02495336 1 0 0
## 31 3.5553481 2.944439 -1.1270118 -0.91130319 -0.94674994 1 0 0
## 32 2.4849066 4.430817 -1.6873995 -1.07002483 0.87546874 2 0 0
## 33 2.4849066 5.288267 -1.4740333 -0.77219039 0.93137637 2 0 0
## 34 2.8903718 5.159055 -1.6044504 -0.72154666 1.11612471 2 0 0
## 35 3.4965076 3.218876 -1.2241755 -0.99155322 -0.65200524 1 0 0
## 36 3.7612001 3.610918 -1.4396951 -0.71539279 -0.73188801 1 0 0
## 37 4.0775374 4.779123 -0.9162907 -0.84863208 1.04133622 1 0 0
## 38 1.6094379 4.204693 -1.7837913 -2.74887220 0.67854103 1 0 0
## 39 3.2958369 4.584967 -1.2765435 -1.11169753 1.18387213 1 0 0
## 40 3.4011974 2.890372 -3.1235656 -1.42711636 -0.98082925 2 0 0
## 41 2.9957323 3.044522 -2.1202635 -1.08175517 -1.19072758 2 0 0
## 42 3.5263605 5.552960 -1.1940225 0.06109510 1.81531322 1 0 0
## 43 3.1780538 4.955827 -1.2909842 -1.21739582 2.39087066 2 0 0
## 44 2.9957323 5.231109 -2.3025851 -0.50087529 1.18814825 3 0 0
## 45 2.0794415 4.779123 -1.4396951 -1.58474530 0.95049890 1 0 0
## 46 0.0000000 4.189655 -0.8627500 -4.71053070 0.09531018 2 0 0
## 47 2.8332133 4.804021 -1.2275827 -0.97286108 0.97531407 1 0 0
## 48 3.4339872 4.890349 -1.2140231 -0.49758040 0.28517894 1 0 0
## 49 2.5649494 3.806662 -0.7787051 -0.76787073 -0.06400533 1 0 0
## 50 3.6375862 5.817111 -0.9675840 -0.40496523 1.28453845 1 0 0
## 51 3.6888795 3.784190 -0.9014021 -0.88673193 -0.43540898 1 0 0
## 52 3.7376696 3.931826 -1.0788097 -0.51082562 -0.01918282 1 0 0
## 53 3.0910425 4.442651 -1.2801342 -0.57270103 -0.05762911 1 0 0
## 54 3.7135721 3.091042 -1.9661129 -0.02429269 0.52591126 1 0 0
## 55 2.7080502 2.302585 -1.7897615 -1.62455155 -2.22562405 1 0 0
## 56 3.2580965 4.343805 -1.2587810 -1.36649173 0.68006194 1 0 0
## 57 3.2580965 4.564348 -1.5511690 -0.69314718 0.22074067 1 0 0
## 58 2.3978953 2.944439 -2.2633644 -2.10373423 -0.35097692 2 0 0
## 59 3.3322045 4.262680 -1.4524342 -1.04982212 0.23744086 1 0 0
## 60 3.5263605 4.477337 -1.3704210 -0.92886951 1.29937389 1 0 0
## 61 3.6888795 3.218876 -1.3356012 -0.91629073 -0.56563386 1 0 0
## 62 3.6635616 2.833213 -2.5257286 -0.44785082 -1.19072758 2 0 0
## 63 3.2958369 3.737670 -1.3547957 3.29583687 -0.23319389 2 0 0
## 64 3.2580965 4.653960 -0.6931472 -0.73188801 0.34995240 1 0 0
## 65 3.5835189 3.401197 -1.4188176 -0.57447565 -0.49102300 1 0 0
## 66 2.7725887 5.141664 -1.3664917 -1.05555280 1.15278477 1 0 0
## 67 2.9444390 3.044522 -1.6873995 -1.49610923 -1.68200861 2 0 0
## 68 3.6375862 4.025352 -2.1803675 -0.88430769 -0.69314718 2 0 0
## 69 3.5835189 3.496508 -2.3025851 -0.44161055 -0.49265832 1 0 0
## 70 2.8332133 5.690359 -1.3356012 -0.83011304 1.05570479 1 0 0
## 71 3.2580965 3.091042 -1.9449106 -1.31676830 -0.34249031 1 0 0
## 72 3.3672958 2.890372 -3.0791139 -1.16475209 -0.77219039 3 0 0
## 73 3.3322045 3.713572 -1.2073117 -1.46967597 0.18730910 1 0 0
## 74 1.6094379 3.044522 -1.2517635 -2.65926004 -1.17441400 1 0 0
## 75 1.7917595 1.098612 -2.9187712 -3.10109279 -2.33304430 5 0 0
## 76 3.9120230 3.912023 -1.9661129 -0.45728486 -0.43232256 1 0 0
## mu_intrvl sig_intrvl
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 0 0
## 7 0 0
## 8 0 0
## 9 0 0
## 10 0 0
## 11 0 0
## 12 0 0
## 13 0 0
## 14 0 0
## 15 0 0
## 16 0 0
## 17 0 0
## 18 0 0
## 19 0 0
## 20 0 0
## 21 0 0
## 22 0 0
## 23 0 0
## 24 0 0
## 25 0 0
## 26 0 0
## 27 0 0
## 28 0 0
## 29 0 0
## 30 0 0
## 31 0 0
## 32 0 0
## 33 0 0
## 34 0 0
## 35 0 0
## 36 0 0
## 37 0 0
## 38 0 0
## 39 0 0
## 40 0 0
## 41 0 0
## 42 0 0
## 43 0 0
## 44 0 0
## 45 0 0
## 46 0 0
## 47 0 0
## 48 0 0
## 49 0 0
## 50 0 0
## 51 0 0
## 52 0 0
## 53 0 0
## 54 0 0
## 55 0 0
## 56 0 0
## 57 0 0
## 58 0 0
## 59 0 0
## 60 0 0
## 61 0 0
## 62 0 0
## 63 0 0
## 64 0 0
## 65 0 0
## 66 0 0
## 67 0 0
## 68 0 0
## 69 0 0
## 70 0 0
## 71 0 0
## 72 0 0
## 73 0 0
## 74 0 0
## 75 0 0
## 76 0 0for (i in 1:nrow(lmen3)) {
lmen3$mu_vol[i] <- mean(c(lmen3$log_vol1[i], lmen3$log_vol2[i], lmen3$log_vol3[i]))
lmen3$sig_vol[i] <- sd(c(lmen3$log_vol1[i], lmen3$log_vol2[i], lmen3$log_vol3[i]))
lmen3$mu_intrvl[i] <- mean(c(lmen3$log_int1[i], lmen3$log_int2[i], lmen3$log_int3[i]))
lmen3$sig_intrvl[i] <- sd(c(lmen3$log_int1[i], lmen3$log_int2[i], lmen3$log_int3[i]))
lmen3$mu_flow[i] <- mean(c(lmen3$log_flow1[i], lmen3$log_flow2[i], lmen3$log_flow3[i]))
lmen3$sig_flow[i] <- sd(c(lmen3$log_flow1[i], lmen3$log_flow2[i], lmen3$log_flow3[i]))
}
lvol3_pred <- as.data.frame(predict(men_vol_model, probs=c(0.1,0.9)))
lint3_pred <- as.data.frame(predict(men_intrvl_model, probs=c(0.1,0.9)))
lflow3_pred <- as.data.frame(predict(men_flow_model, probs=c(0.1,0.9)))
lmen3$lvol3_pred <- lvol3_pred$Estimate
lmen3$lint3_pred <- lint3_pred$Estimate
lmen3$lflow3_pred <- lflow3_pred$Estimate
for (i in 1:nrow(lmen3)) {
lmen3$mu_vol_pred[i] <- mean(c(lmen3$log_vol1[i], lmen3$log_vol2[i], lmen3$lvol3_pred[i]))
lmen3$sig_vol_pred[i] <- sd(c(lmen3$log_vol1[i], lmen3$log_vol2[i], lmen3$lvol3_pred[i]))
lmen3$mu_intrvl_pred[i] <- mean(c(lmen3$log_int1[i], lmen3$log_int2[i], lmen3$lint3_pred[i]))
lmen3$sig_intrvl_pred[i] <- sd(c(lmen3$log_int[i], lmen3$log_int2[i], lmen3$lint3_pred[i]))
lmen3$mu_flow_pred[i] <- mean(c(lmen3$log_flow1[i], lmen3$log_flow2[i], lmen3$lflow3_pred[i]))
lmen3$sig_flow_pred[i] <- sd(c(lmen3$log_flow[i], lmen3$log_flow2[i], lmen3$lflow3_pred[i]))
}
vol3_dat_fit <- fitdist(men_all3$vol3, "lnorm")
vol3_pred_fit <- fitdist(exp(lmen3$lvol3_pred), "lnorm")
int3_dat_fit <- fitdist(men_all3$intrvl3, "lnorm")
int3_pred_fit <- fitdist(exp(lmen3$lint3_pred), "lnorm")
flow3_dat_fit <- fitdist(men_all3$flow3, "lnorm")
flow3_pred_fit <- fitdist(exp(lmen3$lflow3_pred), "lnorm")The graphic posterior predictive check of the predicted variables is in the next chunk. First, the 76 data points and 76 predictions are plotted with the 80% two-sided prediction interval. Second, the data and the predictions are compared with violin plots.
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine#ppcheck vol
p1 <- ggplot() +
geom_point(aes(x=seq(1:76), y=vol3), data=men_all3, shape=4) +
geom_line(aes(x=seq(1:76), y= vol3), data=men_all3, alpha=0.25) +
geom_point(aes(x=seq(1:76), y=exp(Estimate)), data=lvol3_pred, color="blue") +
geom_line(aes(x=seq(1:76), y=exp(Estimate)), data=lvol3_pred, color="blue", alpha=0.25) +
geom_ribbon(aes(x=seq(1:76), ymin=exp(Q10), ymax=exp(Q90)), data=lvol3_pred, color="blue", alpha=0.2) +
scale_y_log10("3rd Void Volume")
# ppcheck interval
p2 <- ggplot() +
geom_point(aes(x=seq(1:76), y=intrvl3), data=men_all3, shape=4) +
geom_line(aes(x=seq(1:76), y= intrvl3), data=men_all3, alpha=0.25) +
geom_point(aes(x=seq(1:76), y=exp(Estimate)), data=lint3_pred, color="green") +
geom_line(aes(x=seq(1:76), y=exp(Estimate)), data=lint3_pred, color="green", alpha=0.25) +
geom_ribbon(aes(x=seq(1:76), ymin=exp(Q10), ymax=exp(Q90)), data=lint3_pred, color="green", alpha=0.2) +
scale_y_log10("3rd Void Interval")
# ppcheck flow
p3 <- ggplot() +
geom_point(aes(x=seq(1:76), y=flow3), data=men_all3, shape=4) +
geom_line(aes(x=seq(1:76), y= flow3), data=men_all3, alpha=0.25) +
geom_point(aes(x=seq(1:76), y=exp(Estimate)), data=lflow3_pred, color="yellow") +
geom_line(aes(x=seq(1:76), y=exp(Estimate)), data=lflow3_pred, color="yellow", alpha=0.25) +
geom_ribbon(aes(x=seq(1:76), ymin=exp(Q10), ymax=exp(Q90)), data=lflow3_pred, color="yellow", alpha=0.2) +
scale_y_log10("3rd Urine Flow Rate")
grid.arrange(p1,p2,p3, nrow=3, ncol=1 )
# Third vol comparison
v1 <- ggplot() +
geom_violin(aes(x=1, y=vol3), data=men_all3) +
geom_label(aes(x=0.5, y=100,
label= paste("lognormal\nmu = ",sprintf("%2.4f",vol3_dat_fit$estimate[1]),
" sigma = ", sprintf("%2.4f", vol3_dat_fit$estimate[2]))),
data=men_all3) +
geom_violin(aes(x=2, y=exp(Estimate)), data=lvol3_pred) +
geom_label(aes(x=2.5, y=30,
label=paste("lognormal\nmu = ", sprintf("%2.4f",vol3_pred_fit$estimate[1]),
" sigma = ",sprintf("%2.4f",vol3_pred_fit$estimate[2]))),
data=lvol3_pred) +
scale_x_continuous("", limits=c(0,3), breaks=c(1,2),
labels=c("NHANES 3rd\nVoid Volume", "Predicted 3rd\nVoid Volume")) +
scale_y_log10("Third Void Volume (ml)")
# thrid interval cmparison
v2 <- ggplot() +
geom_violin(aes(x=1, y=intrvl3), data=men_all3) +
geom_label(aes(x=0.5, y=30,
label= paste("lognormal\nmu = ",sprintf("%2.4f",int3_dat_fit$estimate[1]),
" sigma = ", sprintf("%2.4f", int3_dat_fit$estimate[2]))),
data=men_all3) +
geom_violin(aes(x=2, y=exp(Estimate)), data=lint3_pred) +
geom_label(aes(x=2.5, y=10,
label= paste("lognormal\nmu = ",sprintf("%2.4f",int3_pred_fit$estimate[1]),
" sigma = ", sprintf("%2.4f", int3_pred_fit$estimate[2]))),
data=men_all3) +
scale_y_log10("Third Void Interval (min)") +
scale_x_continuous("", limits=c(0,3), breaks=c(1,2), labels=c("NHANES 3rd\nVoid Interval", "Predicted 3rd\nVoid Interval"))
# third flow comparison
v3 <- ggplot() +
geom_violin(aes(x=1, y=flow3), data=men_all3) +
geom_label(aes(x=0.5, y=1,
label= paste("lognormal\nmu = ",sprintf("%2.4f",flow3_dat_fit$estimate[1]),
" sigma = ", sprintf("%2.4f", flow3_dat_fit$estimate[2]))),
data=men_all3) +
geom_violin(aes(x=2, y=exp(Estimate)), data=lflow3_pred) +
geom_label(aes(x=2.5, y=0.1,
label= paste("lognormal\nmu = ",sprintf("%2.4f",flow3_pred_fit$estimate[1]),
" sigma = ", sprintf("%2.4f", flow3_pred_fit$estimate[2]))),
data=men_all3) +
scale_y_log10("Third Urine Flow (ml/min)") +
scale_x_continuous("", limits=c(0,3), breaks=c(1,2), labels=c("NHANES 3rd\nFlow Rate", "Predicted 3rd\nFlow Rate"))
print(v1)
print(v2)
print(v3)
What we really want is the within person summary statistics
# log mean vol comparison
library(ggthemes)
mu_vol_fit <- fitdist(lmen3$mu_vol, "norm")
mu_vol_pred_fit <- fitdist(lmen3$mu_vol_pred, "norm")
vol_p1 <- ggplot(lmen3) +
geom_violin(aes(x=1, y=mu_vol)) +
geom_label(aes(x=0.5, y=1.8,
label= paste("normal\nmu = ",sprintf("%2.4f",mu_vol_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_vol_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=mu_vol_pred)) +
geom_label(aes(x=2.5, y=1.6,
label= paste("normal\nmu = ",sprintf("%2.4f",mu_vol_pred_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_vol_pred_fit$estimate[2]))),
data=lmen3, size=3) +
scale_x_continuous(limits=c(0,3), breaks=c(1,2), labels=c("data", "pred_volume")) +
scale_y_continuous("log mean volume per individual") +
theme_economist_white()
# sigma vol comparison
sig_vol_fit <- fitdist(lmen3$sig_vol, "gamma")
sig_vol_pred_fit <- fitdist(lmen3$sig_vol_pred, "gamma")
vol_p2 <- ggplot(lmen3) +
geom_violin(aes(x=1, y=sig_vol)) +
geom_label(aes(x=0.5, y=2,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_vol_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_vol_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=sig_vol_pred)) +
geom_label(aes(x=2.5, y=1.5,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_vol_pred_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_vol_pred_fit$estimate[2]))),
data=lmen3, size=3) +
scale_x_continuous(limits=c(0,3), breaks=c(1,2), labels=c("data", "pred_volume")) +
scale_y_continuous("sigma volume per individual") +
theme_economist_white()
grid.arrange(vol_p1, vol_p2, ncol=2)
# clog mean interval comparison
mu_int_fit <- fitdist(lmen3$mu_intrvl, "norm")
mu_int_pred_fit <- fitdist(lmen3$mu_intrvl_pred, "norm")
int_p1 <- ggplot(lmen3) +
geom_violin(aes(x=1, y=mu_intrvl)) +
geom_label(aes(x=0.5, y=3.5,
label= paste("normal\nmu = ",sprintf("%2.4f",mu_int_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_vol_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=mu_intrvl_pred)) +
geom_label(aes(x=2.5, y=3.6,
label= paste("normal\nmu = ",sprintf("%2.4f",mu_int_pred_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_int_pred_fit$estimate[2]))),
data=lmen3, size=3) +
scale_x_continuous(limits=c(0,3), breaks=c(1,2), labels=c("data", "pred_volume")) +
scale_y_continuous("log mean void interval per individual") +
theme_economist_white()
# sigma interval comparison
sig_int_fit <- fitdist(lmen3$sig_intrvl, "gamma")
sig_int_pred_fit <- fitdist(lmen3$sig_intrvl_pred, "gamma")
int_p2 <- ggplot(lmen3) +
geom_violin(aes(x=1, y=sig_intrvl)) +
geom_label(aes(x=0.5, y=2.5,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_vol_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_vol_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=sig_intrvl_pred)) +
geom_label(aes(x=2.5, y=2.6,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_vol_pred_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_vol_pred_fit$estimate[2]))),
data=lmen3, size=3) +
scale_x_continuous(limits=c(0,3), breaks=c(1,2), labels=c("data", "pred_volume")) +
scale_y_continuous("sigma void interval per individual") +
theme_economist_white()
grid.arrange(int_p1, int_p2, ncol=2)
# log mean flow comparison
mu_flow_fit <- fitdist(lmen3$mu_flow, "norm")
mu_flow_pred_fit <- fitdist(lmen3$mu_flow_pred, "norm")
flow_p1 <- ggplot(lmen3) +
geom_violin(aes(x=1, y=mu_flow)) +
geom_label(aes(x=0.7, y=0.6,
label= paste("lognormal\nmu = ",sprintf("%2.4f",mu_flow_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_vol_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=mu_flow_pred)) +
geom_label(aes(x=2.5, y=0.5,
label= paste("lognormal\nmu = ",sprintf("%2.4f",mu_flow_pred_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_int_pred_fit$estimate[2]))),
data=lmen3, size=3) +
scale_x_continuous(limits=c(0,3), breaks=c(1,2), labels=c("data", "pred_volume")) +
scale_y_continuous("log mean flow per individual") +
theme_economist_white()
# sigma flow comparison
sig_flow_fit <- fitdist(lmen3$sig_flow, "gamma")
sig_flow_pred_fit <- fitdist(lmen3$sig_flow_pred, "gamma")
flow_p2 <- ggplot(lmen3) +
geom_violin(aes(x=1, y=sig_flow)) +
geom_label(aes(x=0.5, y=.5,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_flow_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_flow_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2.2, y=sig_vol_pred)) + geom_label(aes(x=2.5, y=1.5,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_flow_pred_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_flow_pred_fit$estimate[2]))),
data=lmen3, size=3) +
scale_x_continuous(limits=c(0,3), breaks=c(1,2), labels=c("data", "pred_volume")) +
scale_y_continuous("sigma flow per individual") +
theme_economist_white()
grid.arrange(flow_p1, flow_p2, ncol=2)
Posterior Predictions for the larger NHANES sample with two urine voids
The values in this larger data set are predicted to obtain values for the log mean and log sd of the volume and interval
##########################################
# predictions from brm
########################################
men_only2 <- read_csv("men_only2a.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## ID = col_double(),
## age = col_double(),
## wt = col_double(),
## ht = col_double(),
## bmi = col_double(),
## int1 = col_double(),
## int2 = col_double(),
## vol1 = col_double(),
## vol2 = col_double(),
## flow1 = col_double(),
## flow2 = col_double()
## )men_only2$grp <- rep(0,1856)
for (i in 1:length(men_only2$int1)) {
if (men_only2$int1[i] <= 120) { men_only2$grp[i] <- 1}
if (men_only2$int1[i] >120 && men_only2$int1[i] <=240) { men_only2$grp[i] <- 2 }
if (men_only2$int1[i] >240 && men_only2$int1[i] <=360) { men_only2$grp[i] <- 3 }
if (men_only2$int1[i] >360 && men_only2$int1[i] <=480) { men_only2$grp[i] <- 4 }
if (men_only2$int1[i] >480) { men_only2$grp[i] <- 5 }
}
lmen2 <- data.frame(
log_age = log(men_only2$age),
log_wt = log(men_only2$wt),
log_ht = log(men_only2$ht),
log_bmi = log(men_only2$bmi),
log_int1 = log(men_only2$int1),
log_int2 = log(men_only2$int2),
log_vol1 = log(men_only2$vol1),
log_vol2 = log(men_only2$vol2),
log_flow1 = log(men_only2$flow1),
log_flow2 = log(men_only2$flow2),
grp = men_only2$grp)
larger_vol_pred <- as_tibble(predict(men_vol_model, newdata = lmen2, allow_new_levels=TRUE))
larger_int_pred <- as_tibble(predict(men_intrvl_model, newdata = lmen2, allow_new_levels=TRUE))
larger_flow_pred <- as_tibble(predict(men_flow_model, newdata = lmen2, allow_new_levels=TRUE))
lmen_added3 <- mutate(lmen2,
log_vol3 = larger_vol_pred$Estimate,
log_int3 = larger_int_pred$Estimate,
log_flow3 = larger_flow_pred$Estimate )
for (i in 1:length(lmen_added3$log_vol1)) {
lmen_added3$mu_vol[i] <- mean(c(lmen_added3$log_vol1[i], lmen_added3$log_vol2[i], lmen_added3$log_vol3[i]))
lmen_added3$sig_vol[i] <- sd(c(lmen_added3$log_vol1[i], lmen_added3$log_vol2[i], lmen_added3$log_vol3[i]))
lmen_added3$mu_intrvl[i] <- mean(c(lmen_added3$log_int1[i], lmen_added3$log_int2[i], lmen_added3$log_int3[i]))
lmen_added3$sig_intrvl[i] <- sd(c(lmen_added3$log_int1[i], lmen_added3$log_int2[i], lmen_added3$log_int3[i]))
lmen_added3$mu_flow[i] <- mean(c(lmen_added3$log_flow1[i], lmen_added3$log_flow2[i], lmen_added3$log_flow3[i]))
lmen_added3$sig_flow[i] <- sd(c(lmen_added3$log_flow1[i], lmen_added3$log_flow2[i], lmen_added3$log_flow3[i]))
} Here, the predictions for between-person variation in log means and log standard deviations for void volume, void interval and urine production rate are fit to normal distributions for the log means and gamma distributions for the log standard deviations.
## mu volume
mu_vol_bigpred_fit <- fitdist(lmen_added3$mu_vol, "norm")
vol2_p1 <- ggplot() +
geom_violin(aes(x=1, y=mu_vol), data=lmen3) +
geom_label(aes(x=0.5, y=2, label = paste("lognormal\nmu = ",sprintf("%2.4f",mu_vol_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_vol_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=mu_vol_pred), data=lmen3) +
geom_label(aes(x=2, y=2,
label= paste("lognormal\nmu = ",sprintf("%2.4f",mu_vol_pred_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_vol_pred_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=3, y=mu_vol), data=lmen_added3) +
geom_label(aes(x=3.5, y=2,
label= paste("lognormal\nmu = ",sprintf("%2.4f",mu_vol_pred_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_vol_pred_fit$estimate[2]))),
data=lmen3, size=3) +
scale_x_continuous(limits=c(0,4), breaks=c(1,2,3), labels=c("mu_vol_data", "pred_mu", "larger_pred_mu")) +
scale_y_continuous("mu volume per individual")
## sigma volume
sig_vol_bigpred_fit <- fitdist(lmen_added3$sig_vol, "gamma")
vol2_p2 <- ggplot() +
geom_violin(aes(x=1, y=sig_vol), data = lmen3) +
geom_label(aes(x=0.5, y=2,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_vol_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_vol_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=sig_vol_pred), data=lmen3) +
geom_label(aes(x=2, y=2,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_vol_pred_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_vol_pred_fit$estimate[2]))),
data=lmen_added3, size=3) +
geom_violin(aes(x=3, y=sig_vol), data=lmen_added3) +
geom_label(aes(x=3.5, y=2,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_vol_bigpred_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_vol_bigpred_fit$estimate[2]))),
data=lmen_added3, size=3) +
scale_x_continuous(limits=c(0,4), breaks=c(1,2,3), labels=c("sig_vol_data", "pred_sig", "larger_pred_sig")) +
scale_y_continuous("sigma volume per individual")
grid.arrange(vol2_p1, vol2_p2, ncol=2)
## mu interval
mu_int_bigpred_fit <- fitdist(lmen_added3$mu_vol, "norm")
int2_p1 <- ggplot() +
geom_violin(aes(x=1, y=mu_intrvl), data=lmen3) +
geom_label(aes(x=0.5, y=3, label = paste("lognormal\nmu = ",sprintf("%2.4f",mu_int_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_int_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=mu_intrvl_pred), data=lmen3) +
geom_label(aes(x=2, y=3, label = paste("lognormal\nmu = ",sprintf("%2.4f",mu_int_pred_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_int_pred_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=3, y=mu_intrvl), data=lmen_added3) +
geom_label(aes(x=3.5, y=3, label = paste("lognormal\nmu = ",sprintf("%2.4f",mu_int_bigpred_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_int_pred_fit$estimate[2]))),
data=lmen_added3, size=3) +
scale_x_continuous(limits=c(0,4), breaks=c(1,2,3), labels=c("mu_int_data", "pred_mu", "larger_pred_mu")) +
scale_y_continuous("mu intrvl per individual")
## sigma interval
sig_int_bigpred_fit <- fitdist(lmen_added3$sig_intrvl, "gamma")
int2_p2 <- ggplot() +
geom_violin(aes(x=1, y=sig_intrvl), data = lmen3) +
geom_label(aes(x=0.5, y=2.3,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_int_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_int_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=sig_intrvl_pred), data=lmen3) +
geom_label(aes(x=2, y=1.8,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_int_pred_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_int_pred_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=3, y=sig_intrvl), data=lmen_added3) +
geom_label(aes(x=3.5, y=1.3,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_int_bigpred_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_int_bigpred_fit$estimate[2]))),
data=lmen_added3, size=3) +
scale_x_continuous(limits=c(0,4), breaks=c(1,2,3), labels=c("sig_int_data", "pred_sig", "larger_pred_sig")) +
scale_y_continuous("sigma interval per individual")
grid.arrange(int2_p1, int2_p2, ncol=2)
## mu flow
mu_flow_bigpred_fit <- fitdist(lmen_added3$mu_flow, "norm")
flow2_p1 <- ggplot() +
geom_violin(aes(x=1, y=mu_flow), data=lmen3) +
geom_label(aes(x=0.5, y=-2.5, label = paste("lognormal\nmu = ",sprintf("%2.4f",mu_flow_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_flow_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=mu_flow_pred), data=lmen3) +
geom_label(aes(x=2, y=-2.5, label = paste("lognormal\nmu = ",sprintf("%2.4f",mu_flow_pred_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_flow_pred_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=3, y=mu_flow), data=lmen_added3) +
geom_label(aes(x=3.5, y=-2.5, label = paste("lognormal\nmu = ",sprintf("%2.4f",mu_flow_bigpred_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",mu_flow_bigpred_fit$estimate[2]))),
data=lmen_added3, size=3) +
scale_x_continuous(limits=c(0,4), breaks=c(1,2,3), labels=c("mu_flow_data", "pred_mu", "larger_pred_mu")) +
scale_y_continuous("mu flow per individual")
## sigma flow
sig_flow_bigpred_fit <- fitdist(lmen_added3$sig_flow, "gamma")
flow2_p2 <- ggplot() +
geom_violin(aes(x=1, y=sig_flow), data = lmen3) +
geom_label(aes(x=0.5, y=1,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_flow_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_flow_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=2, y=sig_flow_pred), data=lmen3) +
geom_label(aes(x=2, y=1,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_flow_pred_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_flow_pred_fit$estimate[2]))),
data=lmen3, size=3) +
geom_violin(aes(x=3, y=sig_flow), data=lmen_added3) +
geom_label(aes(x=3.5, y=1,
label= paste("gamma\nshape = ",sprintf("%2.4f",sig_flow_bigpred_fit$estimate[1]),
"\nrate = ",sprintf("%2.4f",sig_flow_bigpred_fit$estimate[2]))),
data=lmen_added3, size=3) +
scale_x_continuous(limits=c(0,4), breaks=c(1,2,3), labels=c("sig_flow_data", "pred_sig", "larger_pred_sig")) +
scale_y_continuous("sigma flow per individual")
grid.arrange(flow2_p1, flow2_p2, ncol=2)
Here, the smaller sample of 76 men with 3 voids from which the prediction model was developed and the larger sample of men with two voids are combined to create a dataset of 1932 men. This dataset contains males from 6 to 80 years as well as a range o body mass and bmi. From this dataset, selections can be made to obtain a hypothetical population similar to a study population with known biometric characteristics.
men_all <- data.frame(age = c(men_all3$age, men_only2$age),
wt = c(men_all3$wt, men_only2$wt),
ht = c(men_all3$ht, men_only2$ht),
bmi = c(men_all3$bmi, men_only2$bmi),
mu_vol = c(lmen3$mu_vol, lmen_added3$mu_vol),
sig_vol = c(lmen3$sig_vol, lmen_added3$sig_vol),
mu_int = c(lmen3$mu_intrvl, lmen_added3$mu_intrvl),
sig_int = c(lmen3$sig_intrvl, lmen_added3$sig_intrvl),
mu_flow = c(lmen3$mu_flow, lmen_added3$mu_flow),
sig_flow = c(lmen3$sig_flow, lmen_added3$sig_flow)
)
write_csv(men_all, "men_all_020721.csv")
W_mu_vol <- fitdist(men_all$mu_vol, "norm")
summary(W_mu_vol)## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean 3.8093540 0.009004540
## sd 0.3957902 0.006366988
## Loglikelihood: -950.6746 AIC: 1905.349 BIC: 1916.482
## Correlation matrix:
## mean sd
## mean 1 0
## sd 0 1ggplot() +
stat_ecdf(aes(mu_vol), data=men_all, geom="point", size=3, shape=4, alpha=0.2) +
stat_ecdf(aes(rnorm(nrow(men_all), W_mu_vol$estimate[1], W_mu_vol$estimate[2])), geom="step", colour="blue", size=0.9) +
labs(title = "within person log mean void volume", x="log mean void volume", y="cumulative frequency")
W_sig_vol <- fitdist(men_all$sig_vol, "gamma")
summary(W_sig_vol)## Fitting of the distribution ' gamma ' by maximum likelihood
## Parameters :
## estimate Std. Error
## shape 3.419945 0.1051106
## rate 4.493335 0.1487532
## Loglikelihood: -823.6231 AIC: 1651.246 BIC: 1662.379
## Correlation matrix:
## shape rate
## shape 1.0000000 0.9283883
## rate 0.9283883 1.0000000ggplot() +
stat_ecdf(aes(sig_vol), data=men_all, geom="point", size=3, shape=4, alpha=0.2) +
stat_ecdf(aes(rgamma(nrow(men_all), W_sig_vol$estimate[1], W_sig_vol$estimate[2])), geom="step", colour="blue", size=0.9) +
labs(title = "within-person sigma void volume", x="sigma void volume", y="cumulative frequency")
W_mu_int <- fitdist(men_all$mu_int, "norm")
summary(W_mu_int)## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean 4.3042542 0.005861938
## sd 0.2576587 0.004144735
## Loglikelihood: -121.3664 AIC: 246.7328 BIC: 257.8655
## Correlation matrix:
## mean sd
## mean 1.000000e+00 -1.294761e-13
## sd -1.294761e-13 1.000000e+00ggplot() +
stat_ecdf(aes(mu_int), data=men_all, geom="point", size=3, shape=4, alpha=0.2) +
stat_ecdf(aes(rnorm(nrow(men_all), W_mu_int$estimate[1], W_mu_int$estimate[2])), geom="step", colour="blue", size=0.9) +
labs(title = "within-person log mean void interval", x="log mean void void interval", y="cumulative frequency")
W_sig_int <- fitdist(men_all$sig_int, "gamma")
summary(W_sig_int)## Fitting of the distribution ' gamma ' by maximum likelihood
## Parameters :
## estimate Std. Error
## shape 5.607287 0.1753075
## rate 9.134002 0.2987475
## Loglikelihood: -13.16441 AIC: 30.32883 BIC: 41.46145
## Correlation matrix:
## shape rate
## shape 1.0000000 0.9558825
## rate 0.9558825 1.0000000ggplot() +
stat_ecdf(aes(sig_int), data=men_all, geom="point", size=3, shape=4, alpha=0.2) +
stat_ecdf(aes(rgamma(nrow(men_all), shape=W_sig_int$estimate[1], rate=W_sig_int$estimate[2])), geom="step", colour="blue", size=0.9) +
labs(title = "within-person sigma void interval", x="sigma void interval", y="cumulative frequency") 
W_mu_flow <- fitdist(men_all$mu_flow, "norm")
summary(W_mu_flow)## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean -0.4863633 0.011076080
## sd 0.4868438 0.007831823
## Loglikelihood: -1350.713 AIC: 2705.425 BIC: 2716.558
## Correlation matrix:
## mean sd
## mean 1 0
## sd 0 1ggplot() +
stat_ecdf(aes(mu_flow), data=men_all, geom="point", size=3, shape=4, alpha=0.2) +
stat_ecdf(aes(rnorm(nrow(men_all), W_mu_flow$estimate[1], W_mu_flow$estimate[2])), geom="step", colour="blue", size=0.9) +
labs(title = "within-person log mean urine production", x="log mean void urine production", y="cumulative frequency")
W_sig_flow <- fitdist(men_all$sig_flow, "gamma")
summary(W_sig_flow)## Fitting of the distribution ' gamma ' by maximum likelihood
## Parameters :
## estimate Std. Error
## shape 2.330229 0.07027057
## rate 3.086750 0.10383195
## Loglikelihood: -1073.909 AIC: 2151.818 BIC: 2162.951
## Correlation matrix:
## shape rate
## shape 1.0000000 0.8964895
## rate 0.8964895 1.0000000ggplot() +
stat_ecdf(aes(sig_flow), data=men_all, geom="point", size=3, shape=4, alpha=0.2) +
stat_ecdf(aes(rgamma(nrow(men_all), W_sig_flow$estimate[1], W_sig_flow$estimate[2])), geom="step", colour="blue", size=0.9) +
labs(title = "within-person sigma urine productionl", x="sigma urine production", y="cumulative frequency")
#men_all_clean <- men_all[-which(is.na(men_all))]
men_cor <- cor(men_all, method="spearman")
print(men_cor, digits=2)## age wt ht bmi mu_vol sig_vol mu_int sig_int mu_flow
## age 1.000 0.5802 0.458 0.5429 -0.199 0.0311 0.162 0.082 -0.2310
## wt 0.580 1.0000 0.721 0.9355 0.026 -0.0287 0.053 0.135 0.0067
## ht 0.458 0.7206 1.000 0.4747 0.043 -0.0474 0.100 -0.023 -0.0279
## bmi 0.543 0.9355 0.475 1.0000 0.016 -0.0071 0.026 0.190 0.0224
## mu_vol -0.199 0.0261 0.043 0.0164 1.000 0.3623 -0.058 0.090 0.8254
## sig_vol 0.031 -0.0287 -0.047 -0.0071 0.362 1.0000 -0.047 0.176 0.3112
## mu_int 0.162 0.0531 0.100 0.0263 -0.058 -0.0470 1.000 0.444 -0.5657
## sig_int 0.082 0.1345 -0.023 0.1903 0.090 0.1756 0.444 1.000 -0.1643
## mu_flow -0.231 0.0067 -0.028 0.0224 0.825 0.3112 -0.566 -0.164 1.0000
## sig_flow -0.431 -0.3029 -0.329 -0.2352 0.078 0.4176 0.255 0.145 -0.0879
## sig_flow
## age -0.431
## wt -0.303
## ht -0.329
## bmi -0.235
## mu_vol 0.078
## sig_vol 0.418
## mu_int 0.255
## sig_int 0.145
## mu_flow -0.088
## sig_flow 1.000lmen3$sig_vol_ratio <- rep(0,76)
lmen3$sig_int_ratio <- rep(0,76)
for (i in 1:76) { lmen3$sig_vol_ratio[i] <- lmen3$sig_vol[i]/lmen3$sig_vol_pred[i]
lmen3$sig_int_ratio[i] <- lmen3$sig_intrvl[i]/lmen3$sig_intrvl_pred[i] }Comparison with Observational Urinary Behavior Studies
The two studies selected are Smolders et al. (2014) and Parsons et al. (2007). The stated purposes of these studies were to further the understanding of the effects of inter- and intra-individual variation human biomonitoring studies and to investigate bladder diary measurements from asymptomatic individuals during both day and night as well as the relationship with age.
library(readxl)
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionParsons_men <- as_tibble(read_xlsx("Urine paper Dec 2020/Parsons-extract.xlsx", sheet="Men"))## New names:
## * `` -> ...1
## * `` -> ...15
## * `` -> ...16
## * `` -> ...17
## * `` -> ...18Smolders <- read_csv("Urine paper Dec 2020/Smolders.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Subject = col_double(),
## DOW = col_character(),
## Code = col_character(),
## Age = col_double(),
## Gender = col_character(),
## HT = col_double(),
## wt = col_double(),
## Time = col_time(format = ""),
## interval_h = col_double(),
## total_time = col_double(),
## vol = col_double(),
## Cr = col_double(),
## SG = col_double()
## )S2 <- dplyr::filter(Smolders, Subject==2, as.numeric(as.POSIXct(Time))>=25200 && as.numeric(as.POSIXct(Time))<=83000)
S2m <- mean(S2$vol)
S4 <- dplyr::filter(Smolders, Subject==4, as.numeric(as.POSIXct(Time))>=25200 && as.numeric(as.POSIXct(Time))<=83000)
S4m <- mean(S4$vol)
S6 <- dplyr::filter(Smolders, Subject==6, as.numeric(as.POSIXct(Time))>=25200 && as.numeric(as.POSIXct(Time))<=83000)
S6m <- mean(S6$vol)
S8 <- dplyr::filter(Smolders, Subject==8, as.numeric(as.POSIXct(Time))>=25200 && as.numeric(as.POSIXct(Time))<=83000)
S8m <- mean(S8$vol)
df1 <- data.frame(x=rep(3,4), vol=c(S2m, S4m, S6m, S8m))
Pv_fit <- fitdist(Parsons_men$MeanVol_day, "lnorm")
Bv_fit <- fitdist(exp(men_all$mu_vol + 0.5*men_all$sig_vol^2), "lnorm")
Sv_fit <- fitdist(df1$vol, "lnorm")
ggplot() +
geom_violin(aes(x=1, y=MeanVol_day), data=Parsons_men) +
geom_violin(aes(x=2, y=exp(mu_vol + 0.5*sig_vol^2)), data= men_all) +
geom_point(aes(x=x, y=vol), data=df1, shape=4, size=3) +
scale_x_continuous(limits=c(0,4), breaks=c(1,2,3), labels=c("Parsons et al.", "NHANES", "Smolders et al.")) +
scale_y_log10("Mean day_time Void Volume") +
geom_label(aes(x=1, y=10, label=paste("lognormal\nmu = ",sprintf("%2.4f",Pv_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",Pv_fit$estimate[2]))), data=Parsons_men) +
geom_label(aes(x=2, y=10, label=paste("lognormal\nmu = ",sprintf("%2.4f",Bv_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",Bv_fit$estimate[2]))), data= men_all) +
geom_label(aes(x=3, y=10, label=paste("lognormal\nmu = ",sprintf("%2.4f",Sv_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",Sv_fit$estimate[2]))), data=df1)
S2i <- mean(S2$interval_h*60)
S4i <- mean(S4$interval_h*60)
S6i <- mean(S6$interval_h*60)
S8i <- mean(S8$interval_h*60)
df2 <- data.frame(x=rep(3,4), int=c(S2i, S4i, S6i, S8i))
Pi_fit <- fitdist(16*60/Parsons_men$F_day, "lnorm")
Bi_fit <- fitdist(exp(men_all$mu_int + 0.5*men_all$sig_int^2), "lnorm")
Si_fit <- fitdist(df2$int, "lnorm")
ggplot() +
geom_violin(aes(x=1, y=12*60/F_day), data=Parsons_men) +
geom_violin(aes(x=2, y=exp(mu_int + 0.5*sig_flow^2)), data= men_all) +
geom_point(aes(x=x, y=int), data=df2, shape=4, size=3) +
scale_x_continuous(limits=c(0,4), breaks=c(1,2,3), labels=c("Parsons et al.", "NHANES", "Smolders et al.")) +
scale_y_log10("Mean daytime Void Interval") +
geom_label(aes(x=1, y=3000, label=paste("lognormal\nmu = ",sprintf("%2.4f",Pi_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",Pi_fit$estimate[2]))), data=Parsons_men) +
geom_label(aes(x=2, y=3000, label=paste("lognormal\nmu = ",sprintf("%2.4f",Bi_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",Bi_fit$estimate[2]))), data= men_all) +
geom_label(aes(x=3, y=3000, label=paste("lognormal\nmu = ",sprintf("%2.4f",Si_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",Si_fit$estimate[2]))), data=df2)
S2f <- mean(S2$vol/(S2$interval_h*60))
S4f <- mean(S4$vol/(S4$interval_h*60))
S6f <- mean(S6$vol/(S6$interval_h*60))
S8f <- mean(S8$vol/(S8$interval_h*60))
df3 <- data.frame(x=rep(3,4), flow=c(S2f, S4f, S6f, S8f))
Pf_fit <- fitdist(Parsons_men$Flow_day, "lnorm")
Bf_fit <- fitdist(exp(men_all$mu_flow + 0.5*men_all$sig_flow^2), "lnorm")
Sf_fit <- fitdist(df3$flow, "lnorm")
ggplot() +
geom_violin(aes(x=1, y=Flow_day), data=Parsons_men) +
geom_violin(aes(x=2, y=exp(mu_flow + 0.5*sig_flow^2)), data= men_all) +
geom_point(aes(x=x, y=flow), data=df3, shape=4, size=3) +
scale_x_continuous(limits=c(0,4), breaks=c(1,2,3), labels=c("Parsons et al.", "NHANES", "Smolders et al.")) +
scale_y_log10("Mean daytime Urine production") +
geom_label(aes(x=1, y=30, label=paste("lognormal\nmu = ",sprintf("%2.4f",Pf_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",Pf_fit$estimate[2]))), data=Parsons_men) +
geom_label(aes(x=2, y=30, label=paste("lognormal\nmu = ",sprintf("%2.4f",Bf_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",Bf_fit$estimate[2]))), data= men_all) +
geom_label(aes(x=3, y=30, label=paste("lognormal\nmu = ",sprintf("%2.4f",Sf_fit$estimate[1]),
"\nsigma = ",sprintf("%2.4f",Sf_fit$estimate[2]))), data=df3)